####--------------------------------------------####
# Title: Semi-parametric Selection Models for Potentially Non-ignorable 
# Attrition in Panel Studies with Refreshment Samples
# Manuscript ID PA-2013-094
# Y.SI, J.REITER AND S.HILLYGUS
####--------------------------------------------####
#----Bayesian semi-parametric selection models---#
#-------joint for (X Y) & W depends on (X Y)-----#
#----------Plots for APYN data analysis-----------------#
####--------------------------------------------####
###Author: YAJUAN SI
###Latest edit date: 03/30/2014
####--------------------------------------------####
library(ggplot2)
####------------coefficients for AN model------------------------####
beta <-read.csv('output/beta.est.csv',header=FALSE)

var<-c("Intercept","Age: 18-29","Age: 30-44","Age: 45-59","Edu: <High Sch",
       "Edu: High Sch","Edu: Some College","Race: White",
       "Male","Income: <30K","Income: 30-50K","Income: 50-75K",
       "Married","Obama: Favorable (1)","Registered (1)",
     "Interest: Great Deal (1)",
       "Interest: Quite a Bit (1)","Interest: Only Some (1)","Thought: A Lot (1)",
       "Vote: Not Likely (1)","Vote: Some Likely (1)",
       "Vote: Very Likely (1)","Obama: Favorable (2)","Registered (2)",
       "Interest: Great Deal (2)",
       "Interest: Quite a Bit (2)","Interest: Only Some (2)","Thought: A Lot (2)",
       "Vote: Not Likely (2)","Vote: Some Likely (2)",
       "Vote: Very Likely (2)")

beta<-data.frame(cbind(var,beta))
names(beta)<-c("var","median","low","up")
beta$var <- as.character(beta$var)
beta$var <- factor(beta$var, levels=unique(beta$var), ordered=TRUE)
postscript("output/plots/beta.eps")
ggplot() + 
  geom_pointrange(data = beta, aes(x = var, y = median,ymin = low, ymax = up),
                  width = 0.4) +
  coord_flip() +
  geom_point(data = beta, aes(x = var, y = median), colour='red') +
  scale_y_continuous("Posterior summaries for coefficients in the attrition model") +
  scale_x_discrete("") +
  theme(
    plot.background = element_blank()
    ,panel.grid.major = element_blank()
    ,panel.grid.minor = element_blank()
    ,panel.border = element_blank()
    ,panel.background = element_blank()
  )+
  theme(axis.line = element_line(color = 'black'))+
  geom_hline(yintercept=0,linetype = "longdash") +
  theme(axis.text.x = element_text(colour="grey20",size=15,angle=0,hjust=.5,vjust=.5,face="plain"),
        axis.text.y = element_text(colour="grey20",size=18,angle=0,hjust=1,vjust=0,face="plain"),  
        axis.title.x = element_text(colour="grey20",size=20,angle=0,hjust=.5,vjust=0,face="plain"),
        axis.title.y = element_text(colour="grey20",size=0,angle=90,hjust=.5,vjust=.5,face="plain"))
dev.off()


####---------------marginal prob of Y2----------------------####
y2 <-read.csv('output/marprob_y_est.csv',header=FALSE)
jt <-read.csv('output/prob__panelandref_est.csv',header=FALSE)

y2var<-c("Obama Favorability","Registration","Interest","Candidate Thought","Vote Very Likely")
ywvar<-c("Obama Favorability|W","Registration|W","Interest|W","Candidate Thought|W","Vote Very Likely|W")
group<-c(rep("W=1",5),rep("W=0",5))
jt<-data.frame(cbind(group,rep(ywvar,2),jt))
names(jt)<-c("group","ywvar","mean","low","up")
jt$group<-as.character(jt$group)
jt$group <- factor(jt$group, levels=unique(jt$group), ordered=TRUE)
jt$ywvar <- as.character(jt$ywvar)
jt$ywvar <- factor(jt$ywvar, levels=unique(jt$ywvar), ordered=TRUE)
y2<-data.frame(cbind(y2var,y2))
names(y2)<-c("y2var","mean","low","up")
y2$y2var <- as.character(y2$y2var)
y2$y2var <- factor(y2$y2var, levels=unique(y2$y2var), ordered=TRUE)
y2_nw<-cbind(rep("All",5),jt[1:5,2:5])
names(y2_nw)<-c("group","ywvar","mean","low","up")
y2_nw$ywvar <- as.character(y2_nw$ywvar)
y2_nw$ywvar <- factor(y2_nw$ywvar, levels=unique(y2_nw$ywvar), ordered=TRUE)
data_cm<-rbind(jt,y2_nw)
data_cm$ywvar <- as.character(data_cm$ywvar)
data_cm$ywvar <- factor(data_cm$ywvar, levels=unique(data_cm$ywvar), 
 ordered=TRUE)
postscript("outputs/plots/cm.eps")
ggplot()+
  geom_point(data = jt, aes(x = ywvar, y = mean,color=group,shape=group)) +
  geom_pointrange(data = jt, aes(x = ywvar[1:5], y = mean[1:5],
                                 ymin = low[1:5], ymax = up[1:5]),
                  width = 0.4, colour='blue',linetype='dashed') +
  coord_flip() +
  geom_pointrange(data = jt, aes(x = ywvar[6:10], y = mean[6:10],
                                 ymin = low[6:10], ymax = up[6:10]),
                  width = 0.4,colour='orangered3') +
  geom_point(data = y2, aes(x = y2var, y = mean, color='All',shape= 'All'))+
  geom_pointrange(data = y2, aes(x = y2var, y = mean,
                                 ymin = low, ymax = up),colour="black",
                  width = 0.4)+  
  scale_colour_manual(values=c("black","orangered3","blue")) +   
  scale_shape_manual(values=c(5,2,1))+ 
  scale_size_manual(values=c(5.5,5.5,5.5))+
  scale_x_discrete("")+
  scale_y_continuous("Estimates of conditonal and marginal probabilities in the combined data") +
  theme(
    legend.position="bottom",
    legend.direction="horizontal",
    legend.text = element_text(size = 16, face = 'bold'),
    legend.title = element_blank()
  )+
  theme(
    plot.background = element_blank()
    ,panel.grid.major = element_blank()
    ,panel.grid.minor = element_blank()
    ,panel.border = element_blank()
    ,panel.background = element_blank()
  )+
  theme(axis.line = element_line(color = 'black'))+
  theme(axis.text.x = element_text(colour="grey20",size=15,angle=0,hjust=.5,vjust=.5,face="plain"),
        axis.text.y = element_text(colour="grey20",size=20,angle=0,hjust=1,vjust=0,face="plain"),  
        axis.title.x = element_text(colour="grey20",size=20,angle=0,hjust=.5,vjust=0,face="plain"),
        axis.title.y = element_text(colour="grey20",size=0,angle=90,hjust=.5,vjust=.5,face="plain"))
dev.off()

####-------------conditional Y2 in subgroups----------------------####
lv <-read.csv('output/lv3_subgroup_est.csv',header=FALSE)

di <-c("Dem, white, female","Dem, white, male","Dem, nonwhite, female",
       "Dem, nonwhite, male","Rep, white, female","Rep, white, male",
       "Rep, nonwhite, female","Rep, nonwhite, male",
       "Ind, white, female","Ind, white, male","Ind, nonwhite, female",
       "Ind, nonwhite, male")

lv<-data.frame(cbind(di,lv))
names(lv)<-c("di","mean_all","low_all","up_all","mean_panel","low_panel",
             "up_panel","mean_com","low_com","up_com")
mean<-as.numeric(c(lv$mean_all,lv$mean_com))
low<-c(lv$low_all,lv$low_com)
up<-c(lv$up_all,lv$up_com)
type<-c(rep("All",12),rep("Non-attriters",12))
lv1<-data.frame(cbind(type,rep(di,2),mean,low,up))
names(lv1)<-c("type","int","mean","low","up")
lv1$mean<-as.numeric(as.character(lv1$mean))
lv1$low<-as.numeric(as.character(lv1$low))
lv1$up<-as.numeric(as.character(lv1$up))
lv1$int <- as.character(lv1$int)
lv1$int <- factor(lv1$int, levels=unique(lv1$int), ordered=TRUE)
postscript("outputs/plots/lv3.eps")
ggplot()+
  geom_point(data = lv1, aes(x = int, y = mean,color=type,shape=type)) +
  scale_colour_manual(values=c("red","blue")) +   
  scale_shape_manual(values=c(1,2))+ 
  geom_pointrange(data = lv1, aes(x = int[1:12], y = mean[1:12],
                                  ymin = low[1:12], ymax = up[1:12]),
                  width = 0.4, colour='red',linetype='dashed') +
  coord_flip() +
  geom_pointrange(data = lv1, aes(x = int[13:24], y = mean[13:24],
                                  ymin = low[13:24], ymax = up[13:24]),
                  width = 0.4,colour='blue') +
  scale_size_manual(values=c(5.5,5.5))+
  scale_y_continuous("Estimates of probability with campaign interest") +
  scale_x_discrete("") +
  theme(
    legend.position="bottom",
    legend.direction="horizontal",
    legend.text = element_text(size = 16, face = 'bold'),
    legend.title = element_blank()
  )+
  theme(
    plot.background = element_blank()
    ,panel.grid.major = element_blank()
    ,panel.grid.minor = element_blank()
    ,panel.border = element_blank()
    ,panel.background = element_blank()
  )+
  theme(axis.line = element_line(color = 'black'))+
  theme(axis.text.x = element_text(colour="grey20",size=15,angle=0,hjust=.5,vjust=.5,face="plain"),
        axis.text.y = element_text(colour="grey20",size=20,angle=0,hjust=1,vjust=0,face="plain"),  
        axis.title.x = element_text(colour="grey20",size=20,angle=0,hjust=.5,vjust=0,face="plain"),
        axis.title.y = element_text(colour="grey20",size=0,angle=90,hjust=.5,vjust=.5,face="plain"))
dev.off()
